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ABSTRACT 

In  a  wedge  shaped  ocean,  the  method  ot  images  is  used  to 
develop  an  analytical  approximation  of  the  acoustic  pressure 
field.   Contemporary  work  develops  acoustic  doublets  from  a 
combination  of  the  source  and  surface  reflection  image  using 
simple  dipole  theory.   The  method  of  images  is  then  used  to 
sum  the  dipole  images.   This  thesis  matches  dipole  pairs  to 
achieve  a  quadruplet  expansion.   A  computer  program  using  the 
derived  quadruplet  equation  is  then  created  to  verity  the 
results  by  comparing  them  with  the  "URTEXT"  program. 
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I .    INTRODUCTION 

This  thesis  is  a  continuation  in  the  examination  or  the 
acoustic  pressure  field  in  a  wedge  snaped  ocean.   As  the 
focus  of  naval  operations  shifts  to  the  littoral  regions, 
more  attention  is  being  paid  to  the  problem  ot  A.S.W.  in 
acoustically  shallow  water,  i.e.  regions  where  the  sound 
paths  have  multiple  interactions  with  the  surface  ana  bottom. 
To  date,  the  simplest  propagation  model  has  been  the  method 
of  images . 

In  its  present  form,  this  model  assumes  isospeed  water, 
a  pressure  release  upper  surface,  and  a  lossy,  penetrable 
bottom.    for  parallel  surtaces ,  the  source  and  its  images 
form  a  vertical  array  with  the  actual  source  at  the  center. 
Furthermore,  it  tne  source  is  near  trie  surrace  it  can  be 
combined  witn  its  negative  (180°  out  ot  phase)  surtace 
reflection  image  to  form  a  dipole  source.   The  column  ot 
images  can  then  be  considered  a  column  ot  doublet  images. 

The  equations  used  in  this  thesis  are  tne  far  rield 
approximations  of  the  acoustic  doublet.   An  unbalanced 
doublet  is  defined  as  having  two  sources  of  nearly  equal 
amplitude  but  opposite  phase  separated  by  a  distance  d.   It 
the  amplitudes  differ  greatly,  or  if  the  phase  difference  is 
not  180°;  then  the  doublet  approximation  is  not  valid. 

When  the  method  of  images  is  used  in  a  wedge  shaped 
ocean,  the  vertical  column  ot  doublet  images  becomes  a  circle 
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centered  on  the  apex.  In  previous  work  the  upper  and  lower 
image  doublets  were  summed  individually  with  great  success. 
This  thesis  attempts  to  group  upper  and  lower  doublets  into 
double  doublets  or  quadruplets. 

This  method  can  only  be  usea  tor  extremely  small  wedge 
angles.   With  larger  angles  the  sound  paths  from  the  upper 
and  lower  doublet  images  will  be  to  dissimilar,  and  tne 
quadruplet  approximation  will  be  invalid.   Since  typical 
slopes  for  most  of  the  worlds  continental  snelves  are  about 
3°,  a  good  analogy  would  be  that  ot  using  a  ship  mounted 
active  sonar  searching  for  a  diesel  submarine   in  the  shelt 
region . 


FIG.    1 
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II.   WAVE  NUMBERS  AND  SCALING 

Sound  propagation  in  a  shallow  channel  with  a  pressure 
release  surface  parallel  to  a  rigid  bottom  has  a  propagation 
wave  number  for  each  mode  defined  by  the  following  equation 

Where  H  is  the  depth  of  the  channel .   For  a  fast  bottom,  8 
is  the  critical  angle  found  from 


sine--^  U) 


An  analogous  reference  angle  6£  ror  a  slow  bottom  is 

tan0  =— 2  (■*) 

c± 

If  the  depth  H  is  replaced  by  the  scaling  distance, 

and  the  bottom  is  tilted  at  angle  (i ;  tnen  eq .  1  for  the 
cutoff  value  ot  the  lowest  mode  with  a  fast  bottom  becomes 


KX=  2sin6tanp  (b) 


and  for  a  slow  bottom 


KX=  2tan0tanp  (b; 


III.   ACOUSTIC  DOUBLET  FORMATION 

The  sound  propagation  for  spherical  spreading  of  a  point 
source  is 

When  two  point  sources  of  opposite  phase  are  combined  to  form 
a  doublet,  equation  7  can  be  expanded, 

Where  the  "+"  and  "-"  subscripts  refer  to  the  upper  ana  lower 
sources  respectively.   These  values  can  be  expanded 

Where  the  delta  values  indicate  tne  incremental  ditterence 
tor  each  source  compared  to  a  tneoretical  point  source 
located  exactly  between  them.   The  pressure  equation  can  tnen 
be  rearranged  to  yield 

1+AA         1_AA 

p=^£ej^c-kz)   [( ^£_)eJ*Ar_( tS-)e   J'*Ar3         HO) 

r  i-AE       i  +  A£ 

I  I 

An  acoustic  source  located  near  the  surface  has  a 
surface  reflection  image  of  equal  amplitude  and  opposite 
phase.   In  the  far  field  the  image  pair  is  treated  as  a 


single  doublet  source.   For  a  balanced  source,  the  pressure 
equation  can  be  approximated  by 


p=AeJ<*t  [e-jk(z-Az)  _e-JJr(r*Ar)  ]  {li) 


As  seen  in  fig.  2,  Ar  is 


Ar=-^sin6 


( 12; 


Ar-riSinysinO 

(13) 

combining  eqs .  11  and  13,  yields  the  pressure 

P=2j  —  sin(irrisinYsin8)  ej(u>t~kI)  (.14  ) 

To   develop    an    equation    tor    the    unbalanced   doublet,    the 
amplitude    relationships    expressed    in    eq .     9   must    be    rewritten 


A  +A  .        A.-A. 

So    that    for   an   unbalanced   doublet 

P^2-^  [sin^sine)  -j'A^-cos  (Jc^sinO)  ]  ei(at-kr)         (16 ) 

r  Ao 


The  distance  from  the  center  ot  the  doublet  to  the 
receiver  can  be  calculated  using  the  law  of  cosines. 


r2=r12+r22-2r1r2cos(2/3p±6)  U7) 

Where  the  ranges  r^  and  r2  are  measured  from  the  apex  to  the 
source,  and  from  the  apex  to  the  receiver  respectively.   The 
angle  6  is  measured  from  the  surface  down  to  the  receiver. 
The  range  from  eacn  successive  doublet  in  a  wedge  of  angle  (3 
is 


rfl=/r12*r22-2rir2cos(2/iP*d)  U&) 


Where  the  angle  in  the  cosine  term  refers  to  the  angle  or 
each  reflected  image;  those  with  a  "+"  are  the  upper  images, 
and  those  with  a  "-"  are  the  lower  images.   This  equation  can 
be  rearranged  to  yieia 


r=(r.-r2)\l+ — -i-i  (1  -cos  (2np±o)  (^) 

The  first  order  Taylor's  series  expansion  can  tnen  be  used  to 
simplify  eq .  19 

ra-(r1-ra){l+-£l^5  [1-cos (2n$±t) ]}         (20) 

Zl     Z2 


-^ 


L. 


FIG.    2 


IV.   QUADRUPLET  EXPANSION 
The  upper  and  lower  image  doublets  can  be  combined  to  torm 
a  quadruplet.   The  complimentary  doublets  are   180°  out  or  phase, 
and  the  pressure  of  the  new  quadruplet  is 

-sin[icr1Ysin(2/3p-6)]e^'rAr}e-'c;C")r'} 

In  this  case,  both  of  the  interior  exponents  involving  Ar  are 
negative.   The  Ar  terms  are  derived  rrom  a  comparison  to  trie 
range  rrom  the  primary  douoiet,  tnereby  ensuring  that  all  pnase 
angles  are  calculated  with  respect  to  a  common  rererence  point. 
To  find  the  Ar  values  used  in  eq .  21,  an  approximation  or  eq .  2U 
gives 


Ar.+  -r.+  -  {r2-zx)  =  — 1— ?  [1-cos  {2n$±6) }  (22) 


The  cosine  term   can  be  expanded  and  approximated  ror  small 
val ues  or  6  . 

cos  (2/3p±6)  -cos  (2/2p)  ±6sin(2ijP)  l^i 

Equation  21  can  now  be  expanded  into  the  unbalanced  quadruplet 
equation. 


Prr 


[cin[kr1ycin[2nfl+b)]c        a 
-sin[A:rxYsin(27}p-5)  ]  e     I%~^ 


jk-p£-lEin'2i$)  (24; 


r,-r, 


r,-z. 


It  is  convenient  to  define 


^kz^  (l-coa{2n$)  ) 
d=kr16sin(2n$) 


ji=. 


r 


(25) 


So  that  eq .  24  becomes 

Pn=JLLeJ^t-kr)e-j* 

{A„[ sin [Jcrrfsin (220+6)  ]  -7-=~cos[irr,Ysin<2i2p+6)  j  ]  e"J>d      (2b) 
-BD[sin[icr1Y3in(2i2p-o)]  -j -^cos  [Jciyysin  {2flp-6)  ]  ]  cM 

The    values    tor    Aa   and    Ab   must    be    determined    trom    tne    reriectaon 
coefficients . 
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V.   REFLECTION  COEFFICIENTS  FOR  A  SLOW  BOTTOM 
The  reflection  coefficients  used  in  this  thesis  are 
derived  from  the  Rayleigh  reflection  coefficient 

„  i?c-sin0t/sin6v 

*"-= —ix,     •     a  <2'> 


Where 


The  angles  are  the  grazing  angle  Qi   and  angle  or  transmission 
6r   Therefore,  6C  is  defined  as 


Using  trigonometric  relationsnips ,  eq .  it    can  oe 
rewritten . 


bc-Jl  -c2cos28,  /  s  in6 , 

*= *- —  / Ml)) 

br+^/1  -c^cos2Q1/sir\BA 

In   addition,    the    term   under    the    radical    can   be    rewritten. 

v'l-c2cos20i=cv/sin2ei-sin2ec  (31) 

Combining   eqs .    30    and    31,    and    rearranging    yielas 


R_  b-jl-  [sinec/sinej  *_  b-yji-i/x2  32) 


i>+v'l-[sinec/sixi6i]2    Wl"lA*2 
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Where  x  is 


For  a  slow  bottom  with  a  very  small  angle,  eq .  32  can  oe 
approximated  by 


2bc  e, 

The  argument  of  tne  exponential  is 

„_    2bC    _        2p2/Pl 


v'l-c2    j  {c2/ cr)  2-i 


^4; 


(3i>) 


So  that  this  approximation  or  tne  Rayleigh  reflection 
coetticient  can  be  written  in  a  more  concise  torm 

/?=e  aGi  ^b> 

Since  each  ray  intersects  tne  bottom  and  surtace  multiple 
times,  a  product  ot  reflection  coetticients  aiong  eacn 
individual  path  is  required.   This  cumulative  coetticient 
wi 1 1  take  the  r orm 

2C  -1 

&„=  IT  Ri*r&n-a)fk]  U/J 

B"l,3, . . . 

and  the  product  yields 

For  very  small  angles,  this  can  be  approximated  by 
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The  angles  of  incidence  are  needed  to  calculate  tne  proauct. 
The  angles  tor  the  bottom  images  are 

Q0-2{n   1)P  (4U) 

Using  the  law  of  sines,  the  angle  at  the  receiver  is  defined 

e=(i1/r0)sin(2^P  +  o)  ( 41 J 

The  angle  of  incidence  of  tne  ray  trom  the  center  or  tne 
doublet  image  to  tne  receiver  intersecting  the  apparent 
bottom  is 

8=2(/i-I)P+6+€  (4^; 

However,   the  upper  and  lower  images  in  eacn  aoubiet  have 
slightly  different  angles  ot  incidence  with  the  apparent 
bottom 

64=  (222-D  p+5+e4  t4j; 

The  upper  and  lower  angles  at  the  receiver  ror  each  aouhiet 


t±=e±  {z\/±0)lCQS  (2iiP+&)  (44) 

The  individual  angles  of  incidence  are 

8±=60+6  (1+  djLo)  cos  (2/2P+6) )  ±  {i\/i-0)  ycos  (2iip)  (4i>) 

If  the  following  variables  are  used  in  the  coefficients  as 
approximations  ot  the  sine  or  the  angles  ot  incidence 
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ba-^coc2n$ 


(46; 


then    eq .    38    becomes 


Stn  =^e  -*vewri* 


(4?) 


For  each  doublet  image  there  is  an  upper  and  lower  image. 
The  upper  doublet  is  detined  by  AQ  with  the  upper  image  A. 
and  the  lower  image  At .   Likewise,  the  lower  aouoiet  bQ  nas 
the  upper  image  B.  and  tne  lower  image  B  . 


(4b; 


A.    and    AA    detined    by    eq .     15    can    be    combined    with    eq .     48 


Ar=«ncosh  (iibnYO )  e  ns**a       Bn=9Lcosh  (/LbrY« )  ena"5a 
AA=8nsinh  (iLbDya)  e"nM°       AB=8tnsiiih  (/^nY«)  en"M 


(4y; 


When  combined,  they  yield 
AA 


:tanh  (nbnya) 


AB 


tanh  (r?bj«) 


AB         -      Br 
The  result  is  the  full  quadruplet  equation 


(5u) 


D'l        Z 

[e"a"^dsin  (kr1ysin(2n$+6) )  -jtanhb  cos  (ler-LYsin  (2n|i+d)  ) 
-e**^dsin(icr1Ysin(2/2p-5) )  -jLanlib  cos  (ki1ynia{2ii^-tt) )  ] 


(51) 
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FIG.    3 
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VI.   KESULTS 

Three  runs  at  ditferent  ranges  were  made  to  compare  tne 
output  ot  tne  quadruplet  model  witn  mat  or  tne  "UK'lbXT" 
program.   For  each  run,  the  source  angle  is  1°,  and  tne  weage 
angle  is  3°.   The  ratio  ot  the  sound  speed  in  seawater  to 
that  of  the  bottom  is  1.01695,  and  the  ratio  ot  tne  density 
of  seawater  to  that  of  the  bottom  is  0.7.   Tnis  is  equivalent 
to  a  slow  bottom  consisting  ot  a  combination  of  clay  and 
ooze;  conditions  very  common  in  continental  shell  regions. 
The  runs  are  tabulated  at  0.30°  increments  ot  receiver  angle 
from  the  surface  downwards.   The  results  are  snown  in  taoies 
1  through  3.   At  a  receiver  angle  of  zero  ail  three  runs 
result  in  an  amplitude  of  zero,  as  they  snould. 

At  the  shortest  range,  Run  3,  the  results  ot  tne 
quadruplet  model  differ  signi r icant 1 y  trom  those  or  "URTEXT'  . 
The  percentage  airrerence  in  amplitude  is  consistently  aoove 
10%.   With  tne  exception  of  tne  3°  receiver  angle,  tne 
average  percentage  difference  in  amplituae  tor  tne  otner  two 
runs  is  well  below  10%.  Tne  tar  tield  aoubier  approximation 
improves  at  longer  ranges,  and  the  approximations  tor  Ar 
values  become  more  accurate.   Also,  the  difference  in  phase 
angle  is  much  greater  at  the  shortest  range;  half  the  phase 
angles  in  Run  3  differ  by  greater  than  */4,  but  in  the  other 
two  runs,  only  the  phase  angle  at  a  3°  receiver  angle  ditters 
by  more  than  that  amount . 

lb 


Each  ot  these  data  sets  is  accompanied  by   sample  plots 
of  the  quadruplet  pressure  amplitude  and  tne  retiection 
coefficient  for  each  quadruplet.   In  each  case,  the  plots 
were  taken  with  a  receiver  angle  of  0.90°.   The  magnitude  ot 
the  complex  quadruplet  pressure  drops  off  very  quickly  with 
the  quadruplet  number.   It  is  negligible  by  the  5th 
quadruplet.   This  is  a  result  ot  the  effect  ot  tne  rapid 
(exponential)  decay  ot  the  Rayieigh  retiection  coelticient 
with  increased  angle  of  incidence  ot  the  higher  number  ot 
quadrupl et . 

In  practical  terms,  this  means  that  the  quadruplet 
expansion  becomes  less  accurate  witn  larger  receiver  and 
wedge  angles.   As  these  angles  increase,  the  oitierence  in 
the  angles  of  incidence  between  tne  upper  doublet  sound  patn 
and  the  lower  doublet  sound  patn  win  increase 
proportionately.   Wnen  these  ditterences  are  raised  to  the 
exponential  power  in  the  retiection  coet t 1 cienr s ,  tne 
differences  in  amplitude  between  upper  and  lower  oouoiets  at 
the  receiver  will  be  quite  large.   however,  tne  tundamentai 
assumption  in  the  quadruplet  expansion  is  that  the  upper  and 
lower  doublet  images  are  equivalent  in  amplitude,  but 
opposite  in  phase.   Therefore,  the  quadruplet  expansion  is 
only  vaiid  tor  very  small  angles  and  long  ranges. 
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TABLE  1 

BETA=3,GAMMA=1,R1=1,R2=100, DENSITY  RATIO= . 9 , SPEED  RATI0=1 . Oi b95 


DELTA 

URTEXT 

QUAD 

f 
A% 

_ 

1 URTEXT 

QUAD 

1  A° 

AMP. 

AMP. 

AMP. 

|  PHASE 

PHASE 

J  PHASE    'J 

0.30 

0.00010 

0.000092 

8.4 

|  17  .  6 

-23.  9 

i  41 . 5 

0 .  60 

0.00018 

0.000175 

2.8 

|  18.8 

-23  .  5 

1  42.1 

0.90 

0.00024 

0.000243 

1.25 

,19.2 

-22.3 

|  41  .5 

1.20 

0.00032 

0.000288 

10.0 

|  18.3 

-20.8 

j  39  .  1 

1.50 

0.00032 

0.000306 

4.  4 

I  17.5 

-18.  6 

|  3  6.1    i! 

1.80 

0.00030 

U .000295 

1.  7 

1  lb.  2 

-15  .  4 

I  Ji.b    || 

|2.10 

0.00029 

0.000255 

12.0 

|  16.0 

-lO  .  b 

|  -tb.6 

2.40 

0.00020 

0  .000190 

5  .0 

1  15.1 

-2  .  1 

1  1  /  .  I          ,| 

1  2.70 

0.00013 

0.000112 

13  .  9 

|  18.2 

19  .  4 

|  1 .  l          || 

|3.00 

0.00002 

0.000079 

295. 0 

|  2b. 9 

8b  .  b 

1  bl  .  b     || 
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1 
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1 
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TABLE  2 
BETA=3,GAMMA=1,R1=1,R2=50, DENSITY  RAT  10= . 9 , SPEED  RATIU= 1 . 01 o^o 
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|  AMP. 

QUAD 
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AMP. 

F 1 
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1 
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PHASE    | 

A° 

PHASE    || 
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TABLE  3 

BETA=3,GAMMA=1 ,Rl=l ,R2=l0, DENSITY  RAT10= . 9 , SPEED  RAT  10= 1 . U±6*5 


DELTA 

7         1 
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j  AMP. 
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AMP. 

A% 
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1 

QUAD     1 
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APPENDIX  A 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%QUADRUPLET  EXPANSION 

%pquad.m 

%MATLAB 

%Quadruplet  expansion  ot  tne  acoustic 

%pressure  field  in  a  wedge  shaped  ocean. 

%Michael  Joyce 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

B=input(  'wedge  angle  in  degrees^  '  )  ; 

G=mput  (' source    angle    from   surface    in    degrees=  '  j  ; 

D=input (' receiver  angle  from  surtace  in  degrees^  I 

rho=input (' water  to  bottom  density  ratio='); 

cc=input (' water  to  bottom  sound  speed  ratio*'); 

rl=input (' range  of  source  from  apex= ' ) ; 

r2=input (' range  ot  receiver  from  apex=  '  J  ; 

Nl=fix(90/B) ; 
B=B*pi/itfu; 
G=G*pi/180; 
D=D*pi/lbO ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%CALCULATE  SCALING  FACTOR 

2? 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

tb=tan(B) ; 
if  cc<l, 

tl=acos ( cc ) ; 
t2=sin(tl) ; 
else 

tl=acos(l/cc) ; 
t2=tan(tl) ; 
end 
kl=2*t2*tb; 
k=pi/kl ; 

%CALCULATE  CONSTANTS  AND  PRIMARY  DOUBLET 
%%%%%%%%%%%%%%%%%*%%%%%%*%%%%%%%% 

al=2*(l/rho)/(sqrt(cc' 2-1) ) ; 

r=sqrt(rl  2+r2  2-2*rl*r2*cos(D) ) ; 

r3=abs(r2-rl); 

r4=rl*r2/r3  2; 

mu=r2/abs(r2-rl) ; 

q=k*rl*G*D; 

pl=(2/r)*sin(q); 

f  =  0; 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%QUADRUPLET  SUMMATION 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  n=l:l:Nl, 
s  ( n )  =  n ; 

th(n)=2*n*B; 
d(n)=k*rl*D*sin(th(n)); 

phi (n)=K*rl*mux ( 1-cos ( tn( n ) ) ) ; 

^REFLECTION  COEFFICIENTS 
%%%%%%%%%*%%%%%%%%%%%%%%% 

ru(n)=sqrt ( rl  2  +  r2  2- 2*ri * r2*cos  (  th( n ) +D )) 

rl(n)=sqrt(rl  2+r2'2-2*rl*r2*cos(th(n)-D)j 

r5(n)=rl/ru(n) ; 

r6(n)=rl/ rl (nj ; 

aiin^l  +  rsin^cosithin;  )  ; 

bl(n)=r6(n)*cos(th^n) ) ; 

rr(n)  =  -Bxaix(n  2 )  ; 

R(n)=-exp(rr(n)) ; 

a(n)=n*al*D*al(n)  ; 

b(n)=n*al*G*bl(n) ; 
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%PRESSURE 

%%%%%%*%*%*%%%%%*%%%%%%%% 

ul (n)=cosh( b(n) ) ; 
u2(n)  =  exp( -  3* phi (n) ) ; 
u3(n)=(-2/r)*R(n)*ul(nj*u2(nj ; 

%%%%%%%%%%%%%%%%%%%%%%%%% 
%UPPER  IMAGE 

%%%%%%%%%%%%%%%%%%%%*%%%% 

vl(n)=sin( th(n;+D) ; 
v2(n)=k*ri*G*vl(n) ; 
v3(n)=sin( v2(n) ) ; 
v4(n)=cos(v2(n ) ) ; 
v5(n)=-a(n)- }*mu*d( n) ; 
vb(n;=exp( v5 (n ) ) ; 
v7(n)=tanh(b(n) ) ; 
v8ln)=v6(n)Mv3(n)-]*vMnjxv4ln)) 

%%%%%%%%%%%%%%%%%%%%%%%%% 
%LOWER  IMAGE 
%%%%%%%%%%%%%%%%%%%%%%%%% 

wl(n)=sin(th(n)-D) ; 

w2(n)  =  k*rl*G*wl(n/)  ; 
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w3(n)=sin(w2^n)); 
w4(n)=cos(w2(nj ) ; 
w5(n)=a(n)+3*mu*d(n) ; 
w6(n)=exp(w5(n) ) ; 
w7(n)=w6(n)*(w3(n)-;)*v7(n)*w4(n)); 

%%%%%%%%%%%%%%%%%%%%*%%%% 

%SUMMATION 

%%%%%%%%%%%%%%%%%%%%%%%%% 

P(n)=u3(n)*(v8(n)-w7(n;  )  ; 
£=l+p(n); 

pr(n)  =  real ip(n)  )  ; 

pinHn)  =  imag(p(n) ) ; 

pz(n)=sqrt(pr(n)  2  +  p  1  m  (  n ;  2 ) ; 

rnun  )  =  -  ±  x  R  ( n )  ; 
end 
piotvs,pz;;gria; 
xlabel ( ' Quadrupl et  Number'); 
ylabel ( 'Magnatude  or  Quadruplet  Pressure') 
title( 'Quadruplet  Pressure'); 
pause ; 

pi ot ( s , rm) ;grid; 

tit  le(.  '  Rayleigh  Reflection  Coet  ricient  )  ; 
xi abel (' Quadruplet  Number'); 
yiaoen  Coerticient'  ); 
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pause; 

ps=pi+t ; 

ph=imag( ps )/ real ( ps ) ; 

phl=atan(ph) ; 

phase=180*phl/pi 

f  =  f  ; 

Pi=pi; 

ps=ps; 

P  =  sqr t ( rea 1 ( ps  ;  2  +  imagipsj  I ) 

end ; 
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